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Abstract. In recent work Hartmann et al [Phys. Rev. Lett. 102, 057202 
(2009)] demonstrated that the classical simulation of the dynamics of open ID 
quantum systems with matrix product algorithms can often be dramatically 
improved by performing time evolution in the Heisenberg picture. For a closed 
system this was exemplified by an exact matrix product operator solution of 
the time-evolved creation operator of a quadratic fermi chain with a matrix 
dimension of just two. In this work we show that this exact solution can be 
significantly generalized to include the case of an open quadratic fermi chain 
subjected to master equation evolution with Lindblad operators that are linear in 
the fermionic operators. Remarkably even in this open system the time-evolution 
of operators continues to be described by matrix product operators with the 
same fixed dimension as that required by the solution of a coherent quadratic 
fermi chain for all times. Through the use of matrix product algorithms the 
dynamical behaviour of operators in this non-equilibrium open quantum system 
can be computed with a cost that is linear in the system size. We present 
some simple numerical examples which highlight how useful this might be for 
the more detailed study of open system dynamics. Given that Heisenberg picture 
simulations have been demonstrated to offer significant accuracy improvements 
for other open systems that are not exactly solvable our work also provides further 
insight into how and why this advantage arises. 
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1. Introduction 

Developing a more detailed understanding of the numerous intriguing phenomena 
displayed by strongly correlated quantum systems is one of the major theoretical 
challenges in physics today. To meet this challenge a formidable arsenal of non- 
perturbative, renormalization and numerical techniques have been devised. The 
success of these approaches has for the most part been routed in situations at 
or close to equilibrium, while comparatively little is known about the physics of 
strongly-correlated systems far-from-equilibrium. Yet non-equilibrium systems are 
both ubiquitous and of significant practical interest in physics. A typical example is 
where a finite sized strongly-correlated quantum system is driven far-from-equilibrium 
by introducing couplings to several different macroscopic reservoirs, forming an open 
quantum system pQ. Under these circumstances both analytical and numerical 
descriptions of the behaviour of the system become highly non-trivial. 

An immediate need to study open quantum systems is given by the inevitable 
decoherence and dissipation present in any realistic experimental realization of a 
strongly correlated quantum system. Numerous examples of such experiments now 
exist ranging from arrays of Josephson junctions [2], ultra-cold atoms in optical 
lattice O 0] , ion traps [3 El [7] and arrays of coupled microcavities [3 IH [10] . Beyond 
this, however, open systems are becoming increasingly relevant in themselves as efforts 
are made to both understand, and also potentially exploit, the beautiful and subtle 
interplay of the coherent many-body dynamics and incoherent quantum processes 
they possess. One example can be found in quantum information processing [11] 
where the suppression of noise is typically considered a prerequisite. Despite this it is 
found that certain dissipative processes can in fact assist in the preparation of highly 
entangled quantum states [12J, [13l [14]. More generally some of the most common 
occurrences of non-equilibrium physics are in transport problems [15] relevant to 
numerous systems including quantum contacts [16], molecular motors [17] , molecular 
junctions [18] and other low dimensional heat conducting quantum systems |19j . In 
addition to revealing a wealth of non-equilibrium phenomena, including non-diffuse 
heat transfer [20] and negative differential conductance [21], the presence of noise 
has been found, contrary to expectations, to enhance transmission efficiency through 
a dissipative quantum network |22[ 123] where it has a beneficial influence thanks 
to its interplay with destructive quantum interference and energy mismatches |24j . 
It is therefore of technological relevance to better understand quantum mechanical 
effects in driven dissipative strongly-correlated systems in order to exploit them in 
achieving more robust and efficient energy transfer in artificial structures |24] and 
nanomaterials |21j . Finally open quantum systems present a virtually unexplored 
landscape of non-equilibrium phases transitions whose properties are likely to differ 
considerably from conventional equilibrium transitions [25j . 

In this work we shall adopt a master equation [1] description of an open quantum 
system. Our attention is focussed on a specific class of open quantum systems, 
described in detail in Sec. [U which are governed by a quadratic spinless fermionic 
Hamiltonian and coupled to baths described by lineaifj Lindblad operators. Only 
very recently this class of open system was solved semi-analytically [26] with this 
solution later providing strong evidence |27] identifying, somewhat unexpectedly, 
a phase transition in the far-from-equilibrium open XY spin chain with boundary 

| Throughout this paper we use on occasion the phrase linear as a shorthand for operators which 
are 1st order in fermionic creation and annihilation operators, as in Eq. JS}. 
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pumping, but no losses otherwise. Despite being a specialized type of system its 
relevance is elevated by the fact that only a very limited number of exactly solvable 
master equation models are known, namely those involving a single particle, harmonic 
oscillator or spin. Indeed the presence of a non-equilibrium transition in this solvable 
model suggests that it may well come to represent a paradigm for such phenomena 
analogous to the Ising model for quantum phase transitions. 

In this work we present an entirely complementary exact solution of this system 
for the Heisenberg picture evolution of commonly required observables by employing 
a matrix product ansatz 28, 29, 3(J] for the operators, often called a matrix product 
operator (MPO). This result is a significant extension of the exact MPO solution 
presented in [31] for the purely coherent limit of the same system. By using this 
approach our work provides, through the formal structure of the MPO solution, 
further physical insight into why this model is exactly solvable. The utility of our 
solution, however, is only truly revealed when it is combined with powerful matrix 
product based numerical methods, of which density matrix renormalization group 
(DMRG) [321 133] . and more recently its quantum information inspired extension to 
time-evolution [34], [35] [36] , are leading examples. These numerical methods enable 
the exact Heisenberg picture MPO solution of the dynamical evolution of operators to 
be computed under very general situations which are not easily accessible otherwise. 
Given the rich properties of the model [57] this in itself may be very useful. However, 
perhaps of even greater importance, it also provides a significant non-trivial example 
of where Heisenberg picture MPO numerics is exact for an open system. It was 
shown [31] recently that in some cases it is much more efficient and accurate to simulate 
open quantum systems in the Heisenberg picture and thus the underlying numerical 
method used here can be readily applied to more general interacting systems. In 
contrast to the equivalent Schrodinger picture MPO numerics [37] [38] where the 
study of entanglement has provided a crucial understanding of its strengths and 
limitation [39] [40] [41] [42] , the merits of the Heisenberg picture numerics for general 
situations is far less clear. Our result provides further evidence in understanding 
when rigourously exact or good approximate Heisenberg picture MPO solutions may 
exist. We note also that promising results have been found very recently in combining 
solutions of the Heisenberg equations of motion with matrix product representations 
of states for bosonic systems [43] . 

The structure of this paper is as follows. In Sec. [2] we describe the master equation 
and system we shall solve exactly and introduce a particular spin chain model that our 
later numerical calculations will focus on. Our solution exploits the MPO formalism 
and so Sec. [3] describes all the necessary details. We then show in Sec. [4] that the 
Heisenberg picture solution of many operators for a closed coherent system possesses 
an MPO form with a finite dimension, a fact which is crucial for the exact solution 
to be numerically accessible. In Sec. [5] we introduce an ancilla construction which 
reproduces the underlying master equation introduced in Sec. [5] A crucial component 
of this construction is the tracing out of ancillae and the effect of this on an MPO is 
described in Sec. [6] We then combine these observations in Sec.[7jto demonstrate that 
an MPO solution with a bounded dimension, identical to that of the purely coherent 
case, exists for the open system considered. The versatility of this result is highlighted 
in Sec. [5] where we numerically determine the MPO solutions for several situations, 
including the approach to stationarity and a sudden quench of the transverse field. 
Finally in Sec. [9] we conclude and comment on future work. 
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2. Model 

In its most general form the physical system considered in the work is a ID system of 
spinless fermions governed by a quadratic Hamiltonian which reads 



where cj is fermionic creation operator for site j . By demanding H s to be Hermitian 
we can choose a to be real symmetric and b to be real antisymmetric matrices. To 
describe an open fermi lattice we adopt the quantum master equation approach leading 
to a Heisenberg picture evolution of an operator 0(t) that is governed by the Lindblad 
master equation of the form [T] (using % = 1 throughout) 

^-=H{0(t)} + £{0(t)}, (2) 

where the Hamiltonian and Lindblad superoperators are time-independent and defined 
as 

H{0(t)} = i[H s ,0(t)], (3) 
£{0{t)} = Y, ( L lO(t)L 7 - \L\L^O(t) - ±0(t)L\L 7 ) , (4) 

7 

respectively. Here L 1 are Lindblad operators specifying the coupling of the system to 
a set of Markovian baths. We place a restriction on the operators L 7 that they are 
linear in the fermionic creation and annihilation operators with the form 

L 7 = Yj (j-v c 1 + hi c i) ' (5) 

3 

where £ 7 j and l 7 j are complex coefficients. A final constraint, which shall been seen 
in Sec. [5] to be essential to our result, is that 0(t) must be an even ordered operator, 
so PO(i)P = 0(t) where P = JljLiU ~~ ^ c \ c j) IS the parity operator. Since parity is 
conserved by Eq. {2j with quadratic H and linear L 7 's we require only that the initial 
operator O(0) is even. 

Very recently this class of open systems was solved semi-analytically [26j by an 
entirely different approach to that which will be described here. In [26] a sophisticated 
method of constructing a Fock space of operators was employed which maps the 
Liouvillian into a form which can be diagonalized by a procedure analogous to the 
famous solution of the XY Hamiltonian 44J. This solution gives access to a range 
of properties of the system including expectation values of observables for the non- 
equilibrium stationary state and excitations, as well as the spectrum of so-called 
rapidities [55] . We shall exploit this solution later in Sec. [5] for testing the approach 
to stationarity in a dynamical setting. 

The fermionic model outlined has considerable freedom in the non-locality of the 
terms in H s and the operators L 7 . We shall consider a concrete example within this 
class of open fermi systems composed of an XY spin chain with boundary pumping, 
as depicted in Fig. [T] As is well known the Jordan- Wigner transformation 

4=\]\ 1 M [ i <;<» 
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Figure 1. A schematic plot of the open spin chain model considered in this 
work relevant for transport problems. The coherent evolution of the spin chain 
is described by an XY type Hamiltonian H which maps to an effective quadratic 
fermionic Hamiltonian. Boundaries of the chain are subject to couplings to 
baths which are described by Lindblad operators L 1 each of which map to linear 
fermionic operators. 



which relates spin ladder operators a to fermionic creation and annihilation operators 
maps the XY spin-chain Hamiltonian 

N , v N 

3 = 1 V 7 3 = 1 

directly to a spinless fermionic Hamiltonian [44] of the type in Eq. ([I]). Here J is the 
strength of the nearest-neighbour spin coupling, 7 is the anisotropy, B is a transverse 
magnetic field, and er" is the a — {x, y, z} Pauli spin operator on the jth spin. The 
boundary pumping is described by the set of Lindblad operators 

L 2 ^^a^, Li = ^a N . 

L R 

where L + '_ are positive coupling constants. This essentially models a system where 
the two ends of the spin chain are coupled to separate thermal and magnetic baths. 
For an uncoupled chain, where J = 0, the ratios of the local bath couplings 
■pL,R^L,R _ CX p(_2i? /Tl,,r) give the temperature of the thermal state that the 
baths drive the boundary spins to. This spin chain setup is not only of importance 
to heat and spin transport problems [15j in ID but also strong numerical evidence 
suggests it possesses a non-equilibrium phase transition as B is varied |27j . Later in 
Sec. |8] we shall present some exact numerical results for the dynamical behaviour of 
this system possible only through the solution that we will now describe. 

3. Matrix product operators 

The framework in which we cast our exact solution of Eq. j2|) is the matrix product 
representation of operators. Given a system composed of N sites each with a local 
d-dimensional Hilbert space spanned by the states | j) we define the tensor-product 
basis states as | j) = | j%) \ j'2) . . . | jjsr) where j = (ji,j2,---, jw) is a vector of physical 
indices. An arbitrary operator O acting on this system can then be expanded in the 
operator basis | j) (k | as O = J^j I j) s I- A matrix product operator (MPO) is 
where the coefficients o^k of this expansion are expressed in the following form [29|,l30j 

o j:k = (L\ A^ kl A^ 2k2 ...A^ NkN \R), (7) 

where A^™*"™ j s a ma trix, of dimension \ x x, for each site n, selected by 
two independent physical indices j n and k n for that site, while (L | and | R) are 



Exact matrix product solutions in the Heisenberg picture ... 6 




X 

Figure 2. The graphical representation of an MPO for an operator O. For each 
lattice site n the physical indices j n and k n , represented as the thick vertical lines 
respectively, select a \ X \ matrix A' n l J ™' ! ". The joined up thin horizontal lines 
then represent the site-ordered multiplication of these matrices, and finally the 
small circles at the ends depict the boundary vectors (L | and | R) closing the 
chain. 



X-dimensional row and column boundary vectors, respectively. Each expansion 
coefficient oj.k is therefore encoded as a particular ordered product of A matrices 
associated to each site which is contracted to a scalar by the fixed boundary vectors. 

Given that there are in general exponentially many coefficients ojk a matrix 
product representation in Eq. ([7]) yields a highly compact description of an operator 
if it requires only a small dimension \- F° r this reason, and others, matrix product 
representations for both states and operators have been applied with considerable 
success in a variety of related numerical methods. The key to their success is that 
many states or operators, for ID systems at least, can be very accurately approximated 
by a matrix product representation of small dimension despite formally requiring a 
much larger intractable dimension to be exact. In contrast to this the MPO solutions 
we shall present require only a bounded dimension for the representation to be exact 
when describing operators evolving according to the open system introduced in Sec. [2] 
This means that by utilizing one of these matrix product methods, namely the time- 
evolving-block-decimation (TEBD) algorithm, we can evaluate the exact solution 
numerically. However, beyond this much is learnt about the nature of the solution 
by examining the structure of the formal MPO solution itself. For this purpose we 
utilize an entirely lower triangular form for all A-matrices, introduced in [45j [46J [47] , 
which permits exact low-dimensional MPO representations for many operators to be 
constructed easily. The key feature of this approach is that the lower triangular form 
is preserved under the standard matrix product manipulations such as direct sum 
or direct product. This means that if an operator Oa has a MPO representation 
with matrices A of dimension \a, and an operator Ob has one with matrices B with 
dimension \b, then the operator Oa + Ob has matrices A©B and OaOb has matrices 
A ® B with a dimension of at most xa + Xb and xaXb, respectively. Thus much of 
the algebraic convenience of simple product operators (i.e. = 0\®02®---®On 
over a system of N sites and is an MPO with a v = 1) can be extended to highly 
non-trivial operators with an MPO dimension greater than unity. 

For our purposes we need only consider the simplest MPO with a general 2x2 
lower-triangular form. For an operator O we assign the following matrices to each site 

p 
q r 



A = 



where p, q and r are d x d matrices representing local operators on a site. Note that 
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in this compact form of A the physical indices j and k are subsumed into physical 
operators p, q and r, while the row and column indices of the A matrix are the internal 
X = 2 dimensional indices of the MPO representation. To compute the full operator 
described by assigning A to every site we note that the standard multiplication of two 
A matrices is equivalent to the tensor product of the physical operators they contain 
as 



A x A = 



p®p 
) p + r ( 





r <S>r 



For a longer string of multiplications A x A 



x 



A this generalizes to yield 



an operator in the bottom left corner which is the sum of all terms of the form 
r^)---<E>r<S)q^)p<S)---<E)p with the location of the q operator in the string translated. 
Finally for a lower-triangular MPO the full operator is extracted via the left and right 
boundary states (L | = (0, 1) and | R) — (1, 0) T , which select the bottom left operator 
"matrix element" from the matrix product. Using appropriate choices of a, b and c 
many useful single-particle operators can be formed, for example J2j &j is formed by 
each site having a matrix |45j 

1 
, ° Z 1 

Notice that an MPO representation is based on a tensor-product structure and 
therefore implicitly assumes commutativity between local operators appearing in the A 
matrices for different lattice sites. The local operators cannot therefore be fermionic 
directly. For products of such operator sums, which we shall consider shortly, this 
means that MPO's always arrange the resulting local operators in lattice site ordering. 



4. Exact MPO solution for a closed system 

Using the lower-triangular MPO formalism we reexpress the finite-dimensional MPO 
solution described in |31j for any fermionic operator governed by purely coherent 
evolution with a quadratic H s . To do this we need only consider an arbitrary local 
sum of creation and annihilation operators Ci — xci + yc\ . The formal solution to the 
equation of motion of this operator has the standard form 

C t (t)=e nt {Ct} = e iH ' t Cte- iB ' t . 

It can be readily shown that the action of Ti. on Ci for a quadratic H s is 

H{Ct} = i[H.,C t ] 

= ix^2 {a^c,- + h jt c\ } + iy ^ {a^cj + b je c 3 j , 

3 3 

and thus Cg is transformed into a sum of linear operators spread across the lattice. 
The linearity of TL implies that its repeated application any integer number of times 
p as W{Ci\ generates only a linear operator. Now since the formal solution of the 
equation of motion can be expanded as 

oo 

Ce(t) = J2z7H P {Ct}, (8) 



Exact matrix product solutions in the Heisenberg picture 



8 



we establish the well known fact that the Heisenberg picture unitary time evolution of 
the operator Ci governed by a quadratic Hamiltonian is closed. The general solution 
can then be written as 



JV 



(9) 



where ctp,j{t) and Pej(t) are time-dependent complex coefficients containing all the 
non-trivial features of the evolution. To recast this solution in MPO form we apply 
an inverse Jordan- Wigner transformation back to the equivalent spin representation 
giving 



C t (t) = 



N 

E 



n« 

\fc=i 



(10) 



The spin operator on the righthand side can be expressed as a simple 2x2 lower- 
triangular MPO, independent on the number of sites N, with site-dependent matrices 



1,- 







ij(t) 

where X n j(t) — a,ij(t)&j~ + f3ej(t)a~. From our earlier discussion the bottom left 
Xij(t) operator inserts the necessary site and time dependent superposition of spin 
raising and lowering operators into the product, while the bottom right <t| operator 
creates the Jordan- Wigner operator string which establishes the appropriate anti- 
commutative behaviour. 

Since the evolution is unitary the solution to Cj(t) and Cj(t) for all sites j 
automatically provides the time evolution for any string of local sums of creation and 
annihilation operators, i.e. (C p C q . . . Ck)(t) — C p {t)C q (t) . . . Ck{t). Two consequences 
of this are that the dynamics of a quadratic Hamiltonian H conserves the order of any 
initial fcrmionic operator and the MPO solution for the operator string is simply the 
direct product of the MPO solution for each constituent C-operator given in Eq. (fTTj) . 
The latter then straightforwardly determines the fixed matrix product dimension \ 
required for the solution of any given operator string so \ — 2™ for an nth order 
operator. For example, a general quadratic operator C p (t)C q (t) has a 4 x 4 MPO 
representation, independent of p and q given by matrices for each site j as 



(11) 



B0 = 



X qj (t) 













X p j (t) 








X pj {t)X q3 {t) X pj (t)a? 



a]X qj (t) 









(12) 



where X p j(t) and X q j(t) are the site-dependent X operators associated to C p {t) and 
C q {t). The Kronecker product of the MPO solutions gives the appropriately enlarged 

R) = (1, 0, 0, 0) T which select the accumulated 



boundary vectors (L \ — (0, 0, 0, 1) and 
operators in the bottom left corner as 



C p {t)C q {t) 



N 

Y,X pj (t)X qj (t) 

3 = 1 



N 



j>i \k—i 



x qo {t) 
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A general feature of such solutions for strings of C operators is that each constituent 
C operator contributes its own X operator to the representation and shows how highly 
constrained the evolution of operators is in the space of operators, a fact which has 
ultimately permitted such a compact representation. 

Common spin-chain observables such as cr|, ajtrj +1 , and crjcr| +1 are contained in 
this class of 4 x 4 MPO's. Long-range correlations like cr^a^ involve quartic fermionic 
operators, independent of p and q and thus require x — 16. However, the behaviour 
of some operators can be very different. For local spin observables such as crj and 
the fermionic representation obtained via an inverse Jordan- Wigner transformation 
acquires a linearly growing order with the site index j due to the string of (1 — 2c* c) 
operators which appear. Such an operator could then require an exponentially growing 
MPO dimension x = % to describe its exact solution. Correlations like a^a^ behave 
similarly with an exponentially growing dimension x = 4' p ~ 9 ' dependent on their 
separation. What we shall now show in the remainder of this paper is that the x 
required for the MPO solution of even ordered fermionic operators evolving according 
to the open system described in Sec. [2] is identical to that of the purely coherent 
system. 

5. Ancilla master equation construction 

Open quantum systems typically arise when the system of interest interacts with 
a large bath or reservoir, often identified as the system's environment. Using 
this approach Lindblad master equations can be rigorously derived using various 
microscopic models of the system-environment interactions under the Born-Markov 
approximation and in the limit of extremely large reservoirs [T]. To prove that an 
exact finite dimensional MPO representation exists for the open systems introduced 
in Sec. [2] we shall instead employ a derivation of a master equation similar to that of 
non-selective continuous measurement [1 . While this construction itself is perhaps less 
physically motivated it has the advantage for our purposes that it yields a Lindblad 
master equation exactly with no additional approximations. 

A non-selective continuous measurement process involves dividing time into small 
intervals of length St with each interval associated to a separate independent ancilla (or 
probe) forming a time-ordered chain. At the beginning of each interval 8t the system 
evolves coherently and interacts with the associated ancilla which is subsequently 
measured at the end of the interval. Depending on the interaction, measurement and 
ancilla initial state this setup represents a general indirect continuous monitoring of 
the system [351 D] • I n the case where the indirect measurement is ideal the evolution 
of the system is frozen by the quantum-Zeno effect. For more general imperfect 
measurements the system evolves according to a master equation with Hermitian 
Lindblad operators. In order to model the linear fermionic Lindblad operators 
introduced in Eq. ([5]) we modify this construction slightly by considering a different 
class of system- ancilla coupling and trace out rather than measure the ancilla at every 
time step. As we shall show below this setup, depicted in Fig. [3l produces in the 
continuous limit an effective evolution of the system that is again described exactly 
by a Markov master equation with the chain of ancillae representing a manifestly 
delta-correlated environment in time. 

The constructions begins by augmenting the system of N sites with a chain of 
t+1 ancilla sites described by the fermionic modes a t with t = 0, 1, • • • , r. Occupation 
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O '+ & O f+& 

o ' d ~y wo ' 

c ywo »-* trace ^ f -5 f 

system ancilla system ancilla 

Figure 3. (a) The ancilla construction can be visualized as a time-ordered chain 
of ancilla systems associated to each interval of time St and otherwise frozen in 
the initial state | 0) until that particular time. At time t — St the ancilla-system 
Hamiltonian switches the interaction on with the appropriate t — St ancilla for a 
duration St. (b) At time t the Hamiltonian switches the interaction to the next 
ancilla. All earlier ancilla never interact with the system again and can be traced 
out. 



states of the system + ancillae are chosen to be defined by the specific mode ordering 
| n, m) = (c\r ■ ■ ■ (4)"" (4)"* • ' ' (<4r° I vac) , (13) 
where n = (m, • • • , nj\r) and m = (mo, • • • , m T ) are binary vectors of occupation 
numbers over the ancillae and system modes, respectively. By placing the ancillae 
modes to the right this choice, in conjunction with the Jordan- Wigner transformation 
defined in Eq. ensures that any system operator has a spin equivalent of the 
forrri O s (t) = O s (t) ® l a , where l a is the identity over the corresponding ancillae 
spins. This enables the tracing of spins to be completely equivalent to the tracing 
of the corresponding fcrmionic mode. The ancilla mode label t is essentially a time 
label denoting at which time interval the full Hamiltonian of the system will involve 
that ancilla mode, as depicted in Fig. [3J In Eq. (TT3"]) we have also ordered the ancillae 
amongst themselves so their time label increasing inwards from the right so tracing can 
proceed iteratively from the boundary. The full Hamiltonian of the system + ancillae 
is composed of two parts; the time-independent system Hamiltonian H s involving only 
system modes, and Hi(t) which is a time-dependent interaction Hamiltonian between 
the system and ancillae modes. The time-dependence of Hi(t) is taken to be piece- wise 
constant over intervals St giving a full Hamiltonian 

H(t)=H s +Hi(t), (14) 

' t=0 

where S is a system operator and 0(t) is the Heaviside function. Notice that the 
interaction between the lattice and ancilla in Eq. (|14|) depends on 5t and is singular 
in the limit <5t — > 0. This is physically required in order for the ancilla to have a finite 
influence on the lattice in the limit of a vanishingly small interaction time [48 . Also 
the ancillae possess a zero self-Hamiltonian so the only dynamics acting upon them 
is that generated by the terms in Hj(t). As a final definition for this construction we 
take the initial time t — state p of the system + ancillae to have all ancillae modes 
unoccupied, but otherwise arbitrary. 



§ For notational clarity we do not distinguish symbolically between a fermionic operator and its 
Jordan- Wigner transformed spin equivalent. It should be clear from the context which is implied. 
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Let us focus on a particular time t = T6t. Between the time t and t + St the 
Hamiltonian H(t) is time independent and only involves the system modes and the 
ancilla mode o^. For this reason we shall, without loss of generality, restrict our 
considerations to these modes onljfjj To make the connection to MPOs transparent 
we perform a Jordan- Wigner transformation and work with a spin representation. The 
initial density matrix at time t = becomes a spin state in which the lattice and all 
ancilla spins are uncorrelated p = p s ® \ ]) (f | T with p s being an arbitrary spin state 
of the system. At time t the initial operator 0(t) for the system + ancilla spin is 
composed of system modes only and so it transforms to 0(t) = O s {t) (g) It, where 
O s (t) is the system operator resulting from earlier evolution. Similarly the relevant 
interaction terms in H (t) transform to the spin operators 



where P = n^Li a j 1S ^ ne s P m equivalent of the parity operator for the system. During 
the time interval St the formal solution for the evolution is 

0(t + St) = e *{H+H t (t)}St 0{t y- t {H + H z (t)}6t^ (15) 

which will in general leave an operator defined over both the lattice and ancilla spins. 
Since (0{t + 5t)) = tr as [pO(t + St)] the effective time-evolved lattice operator is defined 
by tracing out ancilla spin T resulting in partial expectation value O s (t + St) = 
(t \0(t + St)\ |} T . The expectation value (0{t + St)) can then be expressed solely 
in terms of system operators as (0{t + St)) = tr s [p s O s (t + St)]. Expanding Eq. (TT5)) 
to 2nd order gives 

0{t + St) = 0(t) + iSt[H s + Hi(t),0(t)} 
(iSt) 2 

+ ^ri H s + H i(t), [H. + Hi(t), 0(t)}} + ■■■, 

which can be simplified considerably after the partial expectation value is taken 
due to the special choice of interaction Hi(t) and ancilla initial state | |} T . In 
particular (j \Hi(t)\ |) T = signifying that the ancilla has no direct back action on the 
system [Hill]- The surviving 2nd order term involving Hi{t) is [Hi(t), [Hi(t), 0(t)}] = 
H^tfOit) - 2H l (t)0{t)H l (t) + 0(t)Hi(t) 2 which then also simplifies since 

(T|^W 2 |T) =^S, 



St 

(T \Hi{t)0{t)Ht(t)\ T) = j t S^O s {t)TS : 



Using the resulting evolution 



St 2 

O s {t + St) = O s (t) + iSt[H s , O s {t)} -—[H s , [H.,0,(t)]] (16) 
+ ^ (25 t PO s (t)P5 - S^SOsit) - O a (t)S*S) + ■■■, 

|| Ancillae modes T + 1, ■ ■ ■ , r which are yet to interact are spectators in the proceeding calculation 
since neither 0(t) nor H(t) contain any of these modes. The ancillae modes 0, • ■ ■ , T — 1 which have 
previously interacted may be contained in 0(t). The proceeding calculation is the same regardless of 
whether these modes are traced out before or after the considered time interval St. Thus for brevity 
we assume that they have been traced out before in the same fashion as we shall trace out ancilla 
mode T below. 
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an equation of motion is formed by taking the continuous limit as 




An implicit inverse Jordan- Wigner transformation back to spinless fermions can be 
assumed whereupon we see that the construction has yielded a standard Lindblad 
master equation [49J, [50] with Lindblad operators L = PS. The construction can 
be straightforwardly extended to account for multiple Lindblad operators L 7 by 
introducing more ancillae and additional interactions of the form in Eq. (|14p at each 
time interval. It also admits the option of having explicitly time-dependent Lindblad 
operators £ 7 (t). 

The presence of the F operator relating the coupling S to the resulting Lindblad 
operator L has important consequences. Following our requirements outlined in Seed] 
our aim is for this construction to model linear L operators. For even parity operators 
O(0) the P operator plays no role making the coupling S equivalent to the Lindblad 
operator and therefore linear also. For the same choice of linear coupling S odd parity 
operators O(0) instead evolve according to a different Lindblad superoperator C of 
the form 



with a sign flip of the first term signifying that the Lindblad operators L are now 
higher order. Alternatively, for odd parity operators to evolve with linear Lindblad 
operators L the interaction S must instead include the parity and be higher order. 

6. Tracing out ancilla within an MPO 

A crucial step in the master equation construction is the repeated tracing out of an 
initially uncorrelated ancilla. We now detail the consequences this step has on the 
resulting MPO representation of the system operator O s (t) given an MPO of the full 
operator 0(t). Specifically, let us take 0(t) as being represented by MPO matrices 
A^l for each site j of dimension \, and, without loss of generality, take the ancilla to 
be the last site j = N + 1. 

Given an initial density matrix p s ® p a Fig. 2] shows that the MPO representing 
the system operator O s (t), satisfying tr sa [0(t)p s ® p a ] — tr s [O s (t)p s ], can be found 
by contracting in isolation the single site ancilla density matrix p a with the matrices 
A^- N+1 ^ of 0(t). The contribution of the ancilla to the remaining MPO is then reduced 
to a matrix T = ~J2j k A\ N+1 ^ k {p a )jk whose effect is simply to transform the right 
boundary vector as T| R) = | R'). The MPO for O s (t), defined only over system sites 
j = 1, • • • , N, then retains the same set of matrices for those sites, but possesses 
the new right boundary vector | R'). Thus, so as long as the initial density matrix 
between the system and ancilla is uncorrelated, the dimension of the MPO for O s (t) is 
identical to that of 0(t). This conclusion can be readily seen to hold for any number 
of initially uncorrelated ancilla located at the right edge of the total system. 

For clarity let's consider what type of operators arise from using an arbitrary left 
boundary vector with lower-triangular MPO's. Using the earlier formal solution 
for a generic quadratic operator C p (t)C q (t) given in Eq. (|12|) . the introduction of an 
arbitrary right boundary vector | R') — (a, (3, 7, S) T , while keeping (L \ = (0,0,0,1), 



C{0{t)} 



(S f O(t)S + ^S^SOit) + ±o(t)s^s) , 



(17) 
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Figure 4. Given that the full ancilla-system operator 0(t) is described by an 
MPO the computation of the expectation value (0(t)) is found by contracting it 
with the system + ancilla density matrix and performing a trace which contracts 
away in pairs the remaining physical legs (vertical lines). Since the system and 
ancilla are uncorrelated the trace over the ancilla can be performed in isolation 
reducing its MPO matrices A^"*" 1 ' to a single matrix T. The reduced system 
operator O a (t) retains the same MPO matrices AM spanning its sites j = 1, . . . , N 
but possess a new right boundary vector | R') resulting from transforming | R) with 
the matrix T. 



gives a weighted sum of the bottom row of operator "matrix elements" . More precisely 
it yields an operator 

G{t) = 5 1 + j¥C p (t) + f3¥C q (t) + aC p (t)C q (t), (18) 

which now contains, modulo a parity operator F, linear and zeroth order operators 
that are derived from the operators appearing in the quadratic operator string with 
the admixture determined by the components of | R'). In general an arbitrary right 
boundary vector | R') for an nth order string C p {t)C q (t) . . . Ck{t) generates a sum of 
all operators of order less than or equal to n derived from the constituents of this 
string. When the parent string operator is an even parity operator then lower order 
odd parity terms, such as the linear terms in Eq. (|18p . acquire an additional factor 
of P. The opposite occurs when tracing a parent string operator which has an odd 
parity. 
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The appearance of F operators coincides with terms which do not share the parity 
of the parent string C p (t)C q (t) . . . Ck{t) and thus they violate parity conservation. 
These terms, however, are never generated by the ancilla construction introduced, 
since tracing out the ancilla spin in the specific initial state p a — | T) (T I does not 
generate an arbitrary boundary vector \R'). Instead the allowed structure for | R') 
can be readily discerned by again considering the operator 0{t) = C p (t)C q (t) starting 
with the standard | R) . Tracing out the ancilla spin j = N + 1 results in the partial 
expectation value of the matrix B^^l from Eq. 1|12|) with the state | f ) and yields a 
matrix T as 



where the complex number C ~ (T t) is n °t necessarily zero. Absorbing 

this matrix into the boundary gives R') = (1,0,0, () T signifying that only zeroth 
and quadratic terms can appear. Depending on the nature of the system-ancilla 
interaction tracing out a single ancilla spin for a general nth order operator string can 
in principle generate a boundary vector | R') corresponding to a sum of operators of 
order n, n—2, n—A, . . . , mod(n, 2). This is then consistent with the resulting incoherent 
evolution preserving the parity of the initial fermionic operators. 

This shows that a given lower-triangular MPO's already possess the capacity 
to describe a very specific class of mixed order operators simply by varying one of 
the boundary vectors. As described in Sec. [4] the coherent evolution according to a 
quadratic Hamiltonian H of any one C operator in a string is described by a specific 
time-dependent X operator in its MPO representation. This is true regardless of the 
boundary vector, and so the same coherent quadratic evolution for this type of mixed 
order operator is automatically captured by this MPO solution. 

7. Exact open system MPO solution 

7.1. Building an MPO solution 

The results of the preceding sections can be readily combined to demonstrate that the 
bounded dimension of MPOs seen for coherent quadratic Heisenberg picture evolution 
also applies to even parity operators evolving according to the specific open system 
introduced in Sec. [21 As mentioned in Sec.[5]for even parity operators the requirement 
for linear Lindblad operators is met by using a linear coupling operator S between the 
system and ancilla. This, along with a quadratic system Hamiltonian H s , makes the 
full time-dependent system + ancillae Hamiltonian H{t) quadratic. If all the ancillae 
modes are retained, as depicted in Fig. Ela), then the subsequent evolution would 
be entirely coherent and would represent a purification of the open dynamics of the 
system alone. 

Since we have a coherent quadratic evolution, following the discussion in Sec. [4J 
an exact MPO solution of fixed dimension therefore exists for the full operator 0(t). 
To extract the reduced operator O s {t) for the system for any time < t < rSt the 
entire ancillae chain is traced out, of which only those labelled up to t have any 
relevance. Given that the ancillae and system are uncorrelated initially the tracing 
out of the ancillae has no effect on the resulting MPO dimension for O s (t). The 
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\R(f)) 

Figure 5. A completely coherent representation of the incoherent evolution of 
the open system introduced involves keeping track of the full time-ordered chain 
of ancilla. If at some time < t < rSt the reduced system operator is required 
then the entire ancillae chain is traced out leaving behind a time-ordered product 
of transformation matrices T. The effect of this is to leave the matrices A^' 
describing the system sites unchanged from the form they acquired due to the 
coherent evolution with the ancilla up to time t and instead introduces a time- 
dependent right boundary vector \R(t)). 



tracing out of the ancilla sites yields a product of time-ordered T matricctUJ as shown 
in Fig. [Sib). The incoherent effects induced by the ancilla are entirely captured by a 
time-dependent boundary vector | R(t)) = T t ■ ■ • T,5 t To | R). This therefore establishes 
that for any even-ordered initial system operator O s (0) there is an equality of the 
required MPO dimension for its coherent evolution with a quadratic Hamiltonian and 
its incoherent evolution with this special type of open system. Since mixed order 
operators arising from tracing out any single ancilla can also be coherently evolved, 
with no change in their MPO dimension, this conclusion is independent of when the 
tracing is performed. In particular ancilla may be traced out immediately after they 
interact, as done explicitly in Sec. and thus the bounded MPO dimension applies 
to the continuum limit as well. We demonstrate this with some numerical examples 
in Sec. ® 

7. 2. Properties of the MPO solution 

Here we make some additional comments on the MPO solution found. Firstly, the 
existence of an exact solution for this open system is not simply a consequence of 
the closure of the equations of motion for the lowest order as it is for quadratic 
coherent evolution in Sec. [4] The lack of unitarity of the evolution means that 
(C p Cq . . . Ck)(t) 7^ C p (t)C q (t) . . . Ck(t) in general so knowledge of the evolution of 
lower order operators does not furnish us with knowledge of the evolution of higher 
order products. Secondly, as shown in Sec.[6]the formal structure of the MPO solution 
with a time-varying boundary vector with \ R(t)) ^ (1,0, •■■,0) T implies that the 
evolution of an initial nth order operator will involve a special type of mixed order 



^[ For ancilla related to later times which have yet to interact the transformation matrix T = 1. 
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operator composed of equal or lower order operators only. This behaviour for the 
operator order, revealed by the formal MPO solution, is entirely consistent with what 
is seen for other related open systems [T]. 

A classic example of this is provided by a modified version of the damped 
harmonic oscillator. Here the coherent part of the master equation in Eq. ([3]) has 
H = ujb^b where cu is the oscillator frequency and b is its corresponding bosonic 
annihilation operator. We then take the Lindblad superoperator C in Eq. ((4| as 
being described by a single linear Lindblad operator L = ab + (3¥ analogous to 
the fermionic model introduced in Sec. Since any initial operator of the system 
can be expanded as 0(t) — '^Z nm o m yi{t){b' s ) n (b) m with n, m > and coefficients 
o mn (t), we need only consider the effect of the righthand side of Eq. j| on a 
general term within this expansion. The action of the coherent part is simply 
H{(tf) n (b) m } = iuj(n - m){tf) n (b) m , while the Lindblad contribution gives 

£{(6t)«(6)-} = I(|/3| 2 - H 2 )(n + m)(&t)"(&r 
+ |/3| 2 nm(6 t )"- 1 (6) m ~ 1 

- \a* I3m{m - l){tf) n {b) m - 2 

- ±a(3*n(n-l)(tf) n - 2 (b) m . 

Thus whenever \/3\ > the action of (H + C) is to leave the order of a constituent 
term (b^) n (b) m unchanged as (n + m) or reduced by two. The formal solution 
0(t) = exp(H + £)t{O(0)} for an initial operator O(0) of order n (odd or even) 
will in general include contributions only from orders n, n — 2, n — 4, . . . , mod(n, 2). 
An identical type of analysis can be performed for an open bosonic lattice defined by 
annihilation operators bj for each site j and again governed by a quadratic Hamiltonian 

where a and b are real symmetric matrices, along with linear Lindblad operators 

/- X b J + ^ h) ■ (20) 

3 

This readily confirms that the lack of growth of an operators order seen for a single 
oscillator also applies to the fully bosonic version of the model introduced in Sec. [U 

Applying a similar analysis on the fermionic lattice model itself reveals that for 
the specific Lindblad superoperator C defined in Eq. ^ only initially even ordered 
operators display this closure property. In contrast odd ordered operators can be 
shown to acquire a proliferating order under the repeated application of C. When such 
growth in the order occurs the link between operator order and MPO dimension seen 
for the coherent solution in Sec. 2] suggest that the dimension will not be boundecFI. 
Notice that our ancilla construction applied to odd parity operators does not model 
C, but rather C. Incoherent evolution according to C reverses the situation with odd 
parity operators now displaying no growth in their order. Thus the ancilla construction 
presented in Sec. [5] models an incoherent evolution where all operators 0(t) have a 
bounded order, which in turn permits the bounded dimension MPO solution. It is 
only for even parity operators, however, that this evolution corresponds to the precise 
open system defined in Sec. 

+ Numerical evidence following calculations like those to be presented in Sec. [8] confirms this. 
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Figure 6. (a) The Schmidt coefficients A„ for the splitting between sites 25 
and 26 are plotted on a log-log scale for the MPO's of the operators er| 5 , crfcr| 
and after a time Jt = 50 of evolution. The analytical limit for the MPO 

dimension for each operator is \ = 4, 16, 64, respectively and these are depicted 
by the corresponding dashed vertical lines, (b) The evolution of the central 2- 
magnetization (cfs) (solid) and boundary 2-correlation (""f^fo) (dashed) up to a 
time Jt = 50 starting from a spin-polarized initial state. All these calculations are 
performed on an XY chain of length N = 50 with 7 = 0.75, B/J = 1, F^/J = 0.5, 
T\IJ = 0.3, T«/J = 0.5 and Tf/J = 0.7. 



Finally, the MPO solution offers a more efficient representation than the closure 
of the operator order can provide on its own. In particular once the exact \ is used 
for the MPO its description only grows linearly with the number of sites N. This is 
also an improved scaling compared to alterative approaches to this open fermi system 
exploiting fermionic Gaussian states [52] . They display a N 2 scaling and moreover are 
restricted to considering initial states which are of Gaussian form. The Heisenberg 
picture MPO approach used here can compute properties for any initial state which 
can itself be well approximated by a matrix product state. 

8. Numerical examples 

Having shown that there exists a formal MPO solution with a specific fixed dimension 
X for a given system operator we now show that this solution can be determined 
numerically via Heisenberg picture evolution with the TEBD algorithm [3H 133 131] • 
While the formal solution presented has no restrictions regarding the locality of the 
terms in H s and the Lindblad operators L 7 , efficient integration of the equation of 
motion via the TEBD algorithm requires that terms are nearest neighbour. Moreover 
since linear Lindblad operators involving fermionic creation and annihilation operators 
away from the boundaries acquire a many-spin Jordan- Wigner a z string the numerical 
solution is restricted to noise terms on or one site in from the boundary. This means 
that specific open XY spin chain model introduced in Sec. [2] can be solved with this 
numerical method. 

The numerical solution determined by TEBD unmistakably demonstrates the 
existence of the bounded MPO dimension proven above for this open system, just 
as it does for the coherent limit [S3J. When normalized according to the Frobenius 
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norm, where X)jk°j. k = the MPO solution produced by TEBE0 is in canonical 
form [Ml [SI] 

0jk =^A [11 Y [1]nki A [2] T [2]:i2k2 ...A [N] T [N]jNkN a [a,+11 ^>, (21) 

where the Schmidt decomposition |llj of the operator for any contiguous bipartition 
is explicitly contained in the representation [37]. Here T^- n ^ nkn are sets of matrices, 
different to the lower triangular matrices A["] J " fc ™ used earlier but performing the 
same function. The new important addition to this representation are the diagonal 
matrices A'"' for the bulk (i.e. n = 2, . . . , JV — 1) with diagonal elements equal to the 
Schmidt coefficients A„ of a bipartition after site n. The appropriately sized boundary 
vectors are now ^A^ | = (1,0, ... ,0) and | A[ Ar+1 l) = (1,0,... ,0) T , representing the 
single unit Schmidt coefficient before site 1 and after site N. Once in this canonical 
form Schmidt coefficients allow the effective MPO dimension of the operators to be 
identified by counting the number of significant Schmidt coefficients e < X u < 1, where 
e is a some small threshold. 

For an XY chain of length N = 50 (see Fig. [5] for the parameters used) we have 
calculated the evolution of the operators erf 5 , erf erf and crf4 cr 27- I n Fig- Ela) the 
central Schmidt coefficients for the chosen operators after a time Jt = 50 of evolution 
are shown. A clear cut-off in the A^'s is seen where their value drops in excess of 11 
orders of magnitude. This cut-off is robust to time-evolution and the insignificant A„'s 
are numerical noise that may be safely truncated away. The effective MPO dimension 
given by this cut-off coincides with the dimension expected from the formal solution. 
We may therefore rigidly enforce the exact MPO dimension required and given the lack 
of truncation the only error in the time integration comes from the customary Trotter 
expansion used in TEBD. In Fig. [6jb) the resulting time-evolution of the central z- 
magnetization <r| 5 and boundary-boundary z correlation cr^crfo l& shown for an initial 
spin-polarized state J.j ■ • • j). The transient evolution displays plateaus caused by 
the time it takes for the influence of the boundary pumping to propagate across the 
chain. This is better illustrated in Fig. [7J where the evolution of the z-magnetization 
profile of the entire chain is plotted up to a time Jt = 500. Being plotted with a 
logarithmic timescale it is apparent that the majority of the z-magnetization in the 
bulk is eroded rapidly by the dynamics from its initial value. However, in Fig. [DJa) a 
more detailed comparison of the z-magnetization profile at a time Jt = 500 and the 
stationary profile [SB] reveals that for N = 50 spins this time is only sufficient to drive 
the boundary z-magnetization to their stationary values and that the bulk is still far 
from stationary. Tests reveal that a significantly longer evolution time is needed to 
achieve convergence of the bulk z-magnetizations. 

To demonstrate a time-dependent dynamical scenario^ we consider the simplest 
case of an abrupt quench of the transverse field. Specifically we evolve the initial 
spin-polarized state with B/J = 10 for a time Jt = 500, analogous to the previous 
example. Then at the time Jt — 500 the transverse field is switched instantaneously 
to B/J = 1 and the evolution is continued. In Fig. Ob) the evolution of the central 
z-magnetization (erf 5) is shown as a function of time around the quench point. For 
times Jt < 500 we see that there is a slow change in (crfs) an d it is still quite far 

* The numerical MPO solution will not be lower-triangular. Instead it will be a gauge equivalent 
canonical solution which maintains an orthonormal matrix product structure essential for stability 
and convergence of the numerical algorithm. 

J The application of the TEBD method can be readily adapted to deal with time-dependence in the 
Heisenberg picture. 
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1 

Figure 7. The z-magnetization profile for the entire chain over time plotted 
against r = log 10 (l + Jt) up to a time of Jt = 500. For visual convenience the 
stationary t — » oo magnetization profile determined from the exact solution given 
in |26l is plotted at r = 3. See also Fig. [8{a) for a comparison of the stationary 
profile with that attained after a time Jt = 500. The Hamiltonian parameters 
used are identical to those stated in Fig. \E\ 



from its stationary value of (cfs^ — —0.0161. In the same region of time Fig. [S^b) 
shows (dashed line) the evolution of (a| 5 ) from the previous example where B/J = 1 
throughout which is slightly closer to its stationary value of (cr| 5 ) s = —0.0391 but 
displays a similar rate of convergence. For time Jt > 500, after the quench, there is 
initially a rapid change in (cr| 5 ) which after a time of approximately 15/ J then settles 
down with small oscillates around a new value which again differs from the stationary 
value of the new transverse field. Instead this newly acquired z-magnetization is very 
close to the non-stationary value obtained via constant evolution with B / J = 1. This 
shows that even after a comparatively long evolution time the system has retained a 
significant memory of its initial spin polarized state. 

9. Conclusions 

We have presented a detailed study of the MPO description of a specific class of open 
quantum systems governed by a master equation with a quadratic spinless fermionic 
Hamiltonian and linear fermionic Lindblad operators. By mapping this master 
equation to an entirely coherent quadratic evolution involving additional ancillae we 
have shown that the MPO representation for the evolution of operators with even 
parity possesses a finite and fixed dimension. This has revealed the quadratic nature 
of the evolution underlying this class of master equations and our ancilla construction 
gives decisive insight into why it is exactly solvable. The formal structure of the 
MPO representation also indicates how a given initial operator can evolve into a 
specific type of mixed order operator, consistent with behaviour seen in other simpler 
open systems. Exploiting the fixed MPO dimension the TEBD algorithm allows the 
dynamical evolution of operators in this non-equilibrium open quantum system to be 
computed with a cost that is linear in the system size. The dynamical behaviour 
accessible via the MPO solution presented therefore complements the existing exact 
solution for this models stationary states and spectral properties [26]. We have 
exemplified this by computing some examples involving the approach to stationarity 




Figure 8. (a) A comparison of the stationary (o) ^-magnetization profile (<r|) 
with that attained from a spin-polarized initial state (o) after evolving for a time 
Jt = 500. The central z-magnetization (o"2 5 ) is highlighted with a dashed line, 
(b) For a sudden quench of the traverse field from B/J = 10 to B/J = 1 at a time 
Jt = 500 the time evolution of the central z-magnetization (o"f 5 ) is plotted. The 
dashed line shows the evolution of (cf 5 ) up to a time Jt = 500, already displayed 
in (a) and in Fig. [7] for a constant transverse field B/J = 1. Aside from those 
stated all other Hamiltonian parameters are identical to those in Fig. \E\ 



and the response of the z-magnetization to a sudden quench in the transverse field. 

An interesting calculation, beyond the scope of the current work, is to perform a 
dynamical quenching through the non-equilibrium quantum phase transition. Such a 
dynamical calculation appear to be very demanding with the Schrodinger picture |27j . 
The non-equilibrium transition manifests itself as a discontinuous change in the 
(ofcrj) — (of correlations, but not in other local observables such as energy and 
magnetization. From the MPO perspective of this work the behaviour of this transition 
appears to be very reminiscent of the matrix-product type equilibrium quantum phase 
transitions [54] . Computing a dynamical crossing of this non-equilibrium transition 
could help determine the realistic adiabacity requirements for its observation. 

Beyond this our work has provided an important and non-trivial class of open 
systems with an exact Heisenberg picture MPO representation. This may yet aid 
in determining other models where such solutions exist. For instance it remains to 
be seen whether Heisenberg picture simulability is readily related to the integrability 
of the underlying model [53] . For example a finite sized XXZ chain can be made 
integrable with appropriate boundary fields, however it is not clear that an efficient 
representation exists for commonly required local observables like a z . This raises the 
question as to whether finite-sized MPO representations of certain types of operators 
are possible for systems possessing a Bethe-ansatz solution. This is an interesting 
open problem and would reveal if the MPO formalism can aid in evaluating otherwise 
very complicated quantities from these solutions. For the presently studied XY 
model the non-interacting nature of the effective fermi system for both the open and 
closed system appears to be a crucial property permitting simulability, which is more 
constraining than integrability alone. 

Finally the MPO solution introduced may allow a better understanding of the 
trade-off between efficiencies possible by changing pictures. Future work 55J will look 
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at how quickly the accuracy of Heisenberg picture simulations breakdown when they 
are applied to models which are only weakly perturbed from the exact solution pre- 
sented here. In the context of spin chains the most obvious extensions outside the 
exact solution would be additional (xjcrj +1 interaction terms and/or dephasing noise. 
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